#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
#include "Function.h"
#include "Interpolation.h"

class function1: public Function
{
    public:
        double get_value(const double x)override{
            return 1.0/(1.0+x*x);
        };
}func1;
class function2: public Function
{
    public:
        double get_value(const double x)override{
            return 1.0/(1.0+25*x*x);
        };
}func2;

using namespace std;

int main()
{
int pointNum = 1000;//the number of the output point to plot the function
std::vector<double> outputPoint1(pointNum+1), outputPoint2(pointNum+1), outputPoint3(pointNum+1);
std::vector<double> outputValue1(pointNum+1), outputValue2(pointNum+1), outputValue3_d(pointNum+1), outputValue3_v(pointNum+1);

//problem B
    cout << "problem B: " << endl;
    for(int n1 = 2; n1 <= 8; n1 += 2){
        //define the input point
        std::vector<double> inputPoint1(n1+1);
        std::vector<double> inputValue1(n1+1);
        for (int i = 0; i <= n1; i++) {
            inputPoint1[i] = -5+10.0*i/n1;
            inputValue1[i] = func1.get_value(inputPoint1[i]);
        }
        //compute the output point
        NewtonFormula newtonSolver(n1+1, inputPoint1, inputValue1);
        for (int i = 0; i <= pointNum; i++) {
            outputPoint1[i] = -5+10.0*i/pointNum;
            outputValue1[i] = newtonSolver.getInterpolationValue(outputPoint1[i]);
        }
        //output the result to plot
        string title = "./data/problemB_" + to_string(n1) + ".txt";
        ofstream fout1(title);
        for (int i = 0; i <= pointNum; i++) 
            fout1 << outputPoint1[i] << " " << outputValue1[i] << endl;
        fout1.close();
    }    
    cout << "The result has been output to the file." << endl;
    cout << "==================================================" << endl;
    
    

//problem C
    cout << "problem C: " << endl;
    for(int n2 = 5; n2 <= 20; n2 += 5){
        //define the input point (Chebyshev nodes)
        std::vector<double> inputPoint2(n2);
        std::vector<double> inputValue2(n2);
        for (int i = 0; i < n2; i++) {
            inputPoint2[i] = cos((2*i+1)*M_PI/(2*n2));//Chebyshev nodes
            inputValue2[i] = func2.get_value(inputPoint2[i]);
        }
        //compute the output point
        NewtonFormula newtonSolver(n2, inputPoint2, inputValue2);
        for (int i = 0; i <= pointNum; i++) {
            outputPoint2[i] = -1+2.0*i/pointNum;
            outputValue2[i] = newtonSolver.getInterpolationValue(outputPoint2[i]);
        }
        //output the result to plot
        string title = "./data/problemC_" + to_string(n2) + ".txt";
        ofstream fout2(title);
        for (int i = 0; i <= pointNum; i++) 
            fout2 << outputPoint2[i] << " " << outputValue2[i] << endl;
        fout2.close();
    }
    cout << "The result has been output to the file." << endl;
    cout << "==================================================" << endl;



//problem D
    cout << "problem D: " << endl;
    //define the input point
    std::vector<double> inputPoint3={0.0, 3.0, 5.0, 8.0, 13.0};
    std::vector<double> inputValue3_d={0.0, 225.0, 383.0, 623.0, 993.0};
    std::vector<double> inputValue3_v={75.0, 77.0, 80.0, 74.0, 72.0};
    //(a): predict position and speed when t = 10
    double eps = 1e-8;
    HermiteFormula hermiteSolver(5, inputPoint3, inputValue3_d, inputValue3_v);
    double outputValue3_10_d = hermiteSolver.getInterpolationValue(10.0);
    double outputValue3_10_v = (hermiteSolver.getInterpolationValue(10.0+eps) 
                          - hermiteSolver.getInterpolationValue(10.0-eps))/(2*eps);
    cout << "when t = 10, d = " << outputValue3_10_d << endl;
    cout << "when t = 10, v = " << outputValue3_10_v << endl;
    //(b): whether the velocity exceeds 81
    double maxVelocity = 0, maxVelocityPoint = 0;
    for (int i = 0; i <= pointNum; i++) {
        outputPoint3[i] = i*13.0/pointNum;
        outputValue3_d[i] = hermiteSolver.getInterpolationValue(outputPoint3[i]);
        outputValue3_v[i] = (hermiteSolver.getInterpolationValue(outputPoint3[i]+eps) 
                            - hermiteSolver.getInterpolationValue(outputPoint3[i]-eps))/(2*eps);
        if (outputValue3_v[i] > maxVelocity){
            maxVelocity = outputValue3_v[i];
            maxVelocityPoint = outputPoint3[i];
        }
    }
    if (maxVelocity > 81)
        cout << "The velocity exceeds 81: when t = " << maxVelocityPoint << ", v = " << maxVelocity << endl;
    else
        cout << "The velocity does not exceed 81: when t = " << maxVelocityPoint << ", v = " << maxVelocity << endl;
    //output the result to plot
    ofstream fout3("./data/problemD_d.txt");
    ofstream fout4("./data/problemD_v.txt");
    for (int i = 0; i <= pointNum; i++) {
        fout3 << outputPoint3[i] << " " << outputValue3_d[i] << endl;
        fout4 << outputPoint3[i] << " " << outputValue3_v[i] << endl;
    }
    fout3.close();
    fout4.close();
    cout << "The result has been output to the file." << endl;
    cout << "==================================================" << endl;


    
//problem E
    cout << "Problem E: " <<endl;
    //define the input point
    std::vector<double> inputPoint4={0.0, 6.0, 10.0, 13.0, 17.0, 20.0, 28.0};
    std::vector<double> inputValue4_Sp1={6.67, 17.3, 42.7, 37.3, 30.1, 29.3, 28.7};
    std::vector<double> inputValue4_Sp2={6.67, 16.1, 18.9, 15.0, 10.6, 9.44, 8.89};

    //(a): compute the average weight
    std::vector<double> outputPoint4_ave(pointNum+1), outputValue4_ave_Sp1(pointNum+1), outputValue4_ave_Sp2(pointNum+1);
    NewtonFormula newtonSolverE1(7, inputPoint4, inputValue4_Sp1);
    NewtonFormula newtonSolverE2(7, inputPoint4, inputValue4_Sp2);
    for (int i = 0; i <= pointNum; i++) {
        outputPoint4_ave[i] = 28.0*i/pointNum;
        outputValue4_ave_Sp1[i] = newtonSolverE1.getInterpolationValue(outputPoint4_ave[i]);
        outputValue4_ave_Sp2[i] = newtonSolverE2.getInterpolationValue(outputPoint4_ave[i]);
    }
    //output the result to plot
    ofstream fout5("./data/problemE_a1.txt");
    ofstream fout6("./data/problemE_a2.txt");
    for (int i = 0; i <= pointNum; i++) {
        fout5 << outputPoint4_ave[i] << " " << outputValue4_ave_Sp1[i] << endl;
        fout6 << outputPoint4_ave[i] << " " << outputValue4_ave_Sp2[i] << endl;
    }
    fout5.close();
    fout6.close();
    cout << "The result has been output to the file." << endl;
    
    //(b): predict whether the samples will die after 15 days
    double outputValue4_day43_Sp1 = newtonSolverE1.getInterpolationValue(28.0+15.0);
    double outputValue4_day43_Sp2 = newtonSolverE2.getInterpolationValue(28.0+15.0);
    cout << "The predicted weight of the first sample after 15 days: " << outputValue4_day43_Sp1 << endl;
    cout << "The predicted weight of the second sample after 15 days: " << outputValue4_day43_Sp2 << endl;
    std::vector<double> outputPoint4_pre(pointNum+1), outputValue4_pre_Sp1(pointNum+1), outputValue4_pre_Sp2(pointNum+1);

    for (int i = 0; i <= pointNum; i++) {
        outputPoint4_pre[i] = (28.0+15.0)*i/pointNum;
        outputValue4_pre_Sp1[i] = newtonSolverE1.getInterpolationValue(outputPoint4_pre[i]);
        outputValue4_pre_Sp2[i] = newtonSolverE2.getInterpolationValue(outputPoint4_pre[i]);
    }
    //output the result to plot
    ofstream fout7("./data/problemE_b1.txt");
    ofstream fout8("./data/problemE_b2.txt");
    for (int i = 0; i <= pointNum; i++) {
        fout7 << outputPoint4_pre[i] << " " << outputValue4_pre_Sp1[i] << endl;
        fout8 << outputPoint4_pre[i] << " " << outputValue4_pre_Sp2[i] << endl;
    }
    fout7.close();
    fout8.close();

    cout << "The result has been output to the file." << endl;
    cout << "==================================================" << endl;

}
